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Abstract We study the impact of noise on a neural population rate model of up 
and down states. Up and down states are typically observed in neuronal networks 
as a slow oscillation, where the population switches between high and low fir¬ 
ing rates (Sanchez-Vives and McCormick, 2000). A neural population model with 
spike rate adaptation is used to model such slow oscillations, and the timescale 
of adaptation determines the oscillation period. Furthermore, the period depends 
non-monotonically on the background tonic input driving the population, hav¬ 
ing long periods for very weak and very strong stimuli. Using both linearization 
and fast-slow timescale separation methods, we can compute the phase sensitivity 
function of the slow oscillation. We find that the phase response is most strongly 
impacted by perturbations to the adaptation variable. Phase sensitivity functions 
can then be utilized to quantify the impact of noise on oscillating populations. 
Noise alters the period of oscillations by speeding up the rate of transition be¬ 
tween the up and down states. When common noise is presented to two distinct 
populations, their transitions will eventually become entrained to one another 
through stochastic synchrony. 

Keywords neural adaptation • phase sensitivity function • stochastic synchrony 


1 Introduction 

Cortical networks can generating a wide variety oscillatory rhythms with frequen¬ 
cies spanning five orders of magnitude (Buzsaki and Draguhn, 2004). Slow os¬ 
cillatory activity (0.1-lHz) has been observed in vivo during decreased periods 
of alertness, such as slow wave sleep and anesthesia (Steriade et ah, 1993). Fur¬ 
thermore, such activity can be produced in vitro when bathing cortical slices in 
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a medium with typical extracellular ion concentrations (Sanchez-Vives and Mc¬ 
Cormick, 2000). A key feature of these slow oscillations is that they tend to be an 
alternating sequence of two bistable states, referred to as the up and down states. 
Up states in networks are characterized by high levels of firing activity, due to 
depolarization in single cells. Down states in networks typically appear quiescent, 
due to hyperpolarization in single cells. There is strong evidence that up states 
are neural circuit attractors, that emerge due to synaptic feedback (Cossart et ah, 
2003). This suggests up states may be spontaneous remnants of stimulus-induced 
persistent states utilized for working memory (Wang, 2001) and other network 
computations (Major and Tank, 2004). 

Several different cellular and synaptic mechanisms have been suggested to un¬ 
derlie the transitions between up and down states. One possibility is that the 
network is recurrently coupled with excitation, stabilizing both a quiescent and 
active state (Amit and Brunei, 1997; Renart et ah, 2007). Fluctuations due to prob¬ 
abilistic synapses, channel noise, and randomness in network connectivity can then 
lead to spontaneous transitions between the quiescent and active state (Bressloff, 
2010; Litwin-Kumar and Doiron, 2012). Alternatively, switches between low and 
high activity states may arise by some underlying systematic slow process. For in¬ 
stance, it has been shown that competition between recurrent excitation and the 
negative feedback produced by activity-dependent synaptic depression can lead to 
slow oscillations in firing rate whose timescale is set by the depression timescale 
(Bart et ah, 2005; Kilpatrick and Bressloff, 2010). Excitatory-inhibitory networks 
with facilitation can produce slow oscillations, due to the slow facilitation of feed¬ 
back inhibition that terminates the up state, the down state is then rekindled 
due to positive feedback from recurrent excitation (Melamed et ah, 2008). These 
neural mechanisms utilize dynamic changes in the strength of neural architecture. 
However, Compte et al. (2003) proposed that single cell mechanisms can also shift 
network states between up and down states. The up state is maintained by strong 
recurrent excitation balanced by inhibition, and transitions to the down state occur 
due to a slow adaptation current. Once in the down state, the adaptation current is 
inactivated, and excitation reinitiates the up state. A similar mechanism has been 
utilized in models of perceptual rivalry, where dominance switches between two 
mutually inhibiting populations are due to the build up of a rate-based adaptation 
current (Laing and Chow, 2002; Moreno-Bote et ah, 2007). 

In this paper, we utilize a rate-based model of an excitatory network with spike 
rate adaptation to explore the impact that noise perturbations have upon the rel¬ 
ative phase and duration of slow oscillations. We hnd that, as in the spiking model 
studied by Compte et al. (2003), the interplay between recurrent excitation and 
adaptation produces a slow oscillation in the hring rate of the network. In fact, for 
slow timescale adaptation currents, the oscillations evolve as fast switches between 
a low and high activity state, stable fixed points of the adaptation-free system. 
Since the timescale and slow dynamics of the oscillation are set by the adaptation 
current, we mainly focus on the impact of perturbation to the adaptation variable 
in our model. As we will show, perturbations of the activity variable have much 
lower impact on the oscillation phase. Introducing noise into the adaptation vari¬ 
able of the population model leads to a speeding up of the slow oscillation, due to 
early switching between the low and high activity state. 

Another remarkable feature of slow oscillations, observed during slow-wave 
sleep and anesthesia, is that the up and down states tend to be synchronized 
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across different regions of cortex and thalamns (Steriade et al., 1993; Massimini 
et al., 2004). Specifically, both the np and down states start near synchronously 
in cells located up to 12nim apart (Volgushev et ah, 2006). Such remarkable co¬ 
herence between distant network activity cannot be accomplished by single cell 
mechanisms, but require either long range network connectivity or some external 
signal forcing entrainment (Traub et ah, 1996; Smeal et ah, 2010). Activity tends 
to originate from several different foci in the network, quickly spreading across 
the rest of the network on a timescale orders of magnitude faster than the oscil¬ 
lation itself (Compte et ah, 2003; Massimini et ah, 2004). The fact that the onset 
of quiescence is fast and well synchronized means there must be either a rapid 
relay signal between all foci or there is some global signal cueing the down state. 
Rather than suggest a disynaptic relay, using long range excitation acting on local 
inhibition, we suggest that background noise can serve as a synchronizing global 
signal (Ermentrout et ah, 2008). Noisy but correlated inputs have been shown to 
be capable of synchronizing uncoupled populations of phase oscillators (Teramae 
and Tanaka, 2004) as well as experimentally recorded cells in vitro (Galan et ah, 
2006). Here we will show correlated noise is a viable mechanisms for coordinating 
slow oscillations in distinct uncoupled neural populations. 

The paper is organized as follows. We introduce the neural population model 
in section 2, indicating the way external noise is incorporated into the model. 
In section 3, we demonstrate the periodic solutions that emerge in the noise-free 
model, demonstrating it is possible to derive analytical expressions for the oscil¬ 
lation period in the case of steep firing rate functions. Then, in section 4 we show 
how to derive phase sensitivity functions that describe how external perturbations 
to the periodic solution impact the asymptotic phase of the oscillation. As demon¬ 
strated, the impact of perturbations to the adaptation variable is much stronger 
than activity variable perturbations, especially for longer adaptation timescales. 
Thus, our studies of the impact of noise mainly focus on the effects of fluctuations 
in the adaptation variable. We find, in section 5, that adding noise to the adap¬ 
tation variable leads to up and down state durations that are shorter and more 
balanced, so that the up and down state last for similar lengths of time. In sec¬ 
tion 6, we demonstrate that slow oscillations in distinct populations can become 
entrained to one another when both populations are forced by the same common 
noise signal. This phenomenon is robust to the introduction of independent noise 
in each population, as we show in section 7. 


2 Adaptive neural populations: deterministic and stochastic models 

We begin by describing the models we will use to explore the impact of external 
perturbations on slow oscillations. Motivated by Compte et al. (2003), we will focus 
on a neural population model with spike rate adaptation, akin to mutual inhibitory 
models used to study perceptual rivalry (Laing and Chow, 2002; Moreno-Bote 
et ah, 2007). 

Single population model. In a single population, neural activity u{t) receives 
negative feedback due to a subtractive spike rate adaptation term (Benda and 
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Herz, 2003) 


u{t) = —u{t) + f{au{t) — a{t) + 7), (la) 

Ta{t) = —a{t) + 4>u{t). (lb) 


Here, u represents the mean firing rate of the neural population with excitatory 
connection strength a. The negative feedback variable a is spike frequency adap¬ 
tation with strength tp and time constant t. For some of our analysis we will utilize 
the assumption t ^ 1, based on the fact that many form for spike rate adaptation 
tend to be much slower than neural membrane time constants (Benda and Herz, 
2003). The constant tonic drive I initiates the high hring rate (up) state, and slow 
adaptation eventually attenuates activity to a low firing rate (down) state. Weak 
but positive drive 7 > 0 is meant to model the presence of low spiking threshold 
cells that spontaneously hre, utilized as a mechanism for initiating the up state 
in Compte et al. (2003). The hring rate function / is monotone and saturating 
function such as the sigmoid 

= 1 + e-y. ■ (2) 

Commonly, in studies of neural held models, the high gain limit of (2) is taken to 
yield the Heaviside hring rate function (Amari, 1977; Laing and Chow, 2002) 



1 : a; > 0, 
0 : a; < 0, 


( 3 ) 


which often allows for a more straightforward analytical study of model dynamics. 
We exploit this fact extensively in our study. Nonetheless, we have also carried 
out many numerical simulations of the model for a smooth hring rate function 
(2), and they correspond to the results we presentfor sufficiently high gain. Note, 
this form of adaptation is often referred to as subtractive negative feedback, as 
current is subtracted from the population input. Alternative models of slow neural 
population oscillations have employed short term synaptic depression (Tabak et ah, 
2000; Bart et ah, 2005; Kilpatrick and Bressloff, 2010), a form of divisive negative 
feedback. 

A primary concern of this work is the response of (1) to external perturbations, 
acting on the activity u and adaptation a variables. To do so, we will use both 
an exact method and a linearization to identify the phase response curve of the 
limit cycle solutions to (1). Understanding the susceptibility of limit cycles (1) to 
inputs will help us understand ways in which noise will inhuence the frequency 
and regularity of oscillations. 

Stochastic single population model. Following our analysis of the noise- 
free system, we will consider how huctuations inhuence oscillatory solutions to (1). 
To do so, we will employ the following Langevin equation for (1) forced by white 
noise 


du(t) = [-u{t) + f{au{t) - a{t) + 7)] dt + d^„(t) (4a) 

da(t) = [-a(t) -h (j>u{t)] dt/r + d^a{t), (4b) 


where we have introduced the independent Gaussian white noise processes ^u(t) 
and £,a(t) with zero mean {^u{t)) = (Ca(t)) = 0 and variances (Cii(t)^) = Uut 
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and = cTat- Extending our results concerning the phase response curve, 

we will explore how noise forcing impacts the statistics of the resulting stochastic 
oscillations in (4). In particular, since we find noise tends to impact the phase of 
the oscillation more strongly when applied to the adaptation variable, we will tend 
to focus on the case = 0. 

Stochastic dual population model. Finally, we will focus on how correla¬ 
tions in noise-forcing impact the coherence of two distinct uncoupled populations 


dui = [-Ml (t) + f{aui {t) - ai {t) + /)] dt -|- d,*^ (5a) 

dai = [-ai(t)-h Mi(t)] dt/r-I-d^a (5b) 

du2 = [-U2{t) + f{au2{t) - 02(t) -h /)] dt -h d,*u (5c) 

da2 = [-02(1) -h U2{t)] dt/r + d,*a- (5d) 


Thus, the system (5) describes the dynamics of two distinct neural populations 
Ml and M 2 , with inputs I. Our main interest lies in the impact the noise terms 
have upon the phase relationship between the two systems’ states. In this version 
of the model, noise to the activity variables is totally correlated, as is noise 
to the adaptation variables ^a- Thus, all means are zero and (fu(l)) = cTu = 
D rit. Furthermore, (^ct(^)) — — l-^ 22 ^. For this study, we assume there are no 

correlations between activity and adaptation noise, so (Cu(t)Ca(l)) = 0. A more 
general version of the model (5) would consider the possibility of independent noise 
in each population 

dui = [-Mi(t) + f{aui{t) - ai{t) + I)] dt -|- Xud^uc -h \/l - xSd^tii (6a) 

dfll = [-Ml {t) -h Ml (t)] dt/T + Xad^ac + \/1 “ Xadfal (6b) 

dM 2 = [-U 2 {t) -h f{au 2 {t) - 02 ( 1 ) + /)] dt + Xud^uc -h (6c) 

da 2 = [-02 (t) -h U 2 {t)] dt/r + Xadfac + Vl - Xadfa2. (6d) 

Noise terms all have zero mean and variances dehned = Dujt and 

= Dajt {j = l,2,c). To ease calculations, we take D„i = D „2 = 
Dim = cTu and Dai = Da 2 = Dai = Oa- The degree of noise correlation between 
populations is controlled by the parameters Xu and Xa, so in the limit Xu,a —t 1, 
the model (6) becomes (5). 


3 Periodic solutions of a single population 

We begin by studying periodic solutions of the single population system (1), as 
demonstrated in Fig. lA. First, we note that for hring rate functions / with finite 
gain, we can identify the emergence of oscillations by analyzing the stability of the 
equilibria of (1). That is, we assume (m, a) = (0,0), so the system becomes 

M = f{au — a + I), 
a = (j}U, 

which can be reduced to the single equation 


u = f{{a - (l>)u + 1) = g{u). 


( 7 ) 
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Fig. 1 Single adapting neural population (1) generates slow oscillations. (A) Numerical simu¬ 
lation of (1) for adaptation timescale r = 100 (Is) and input I = 0.2. (B) Partitioning of (r, /) 
parameter space shows the range of inputs I leading to oscillations expands as the adaptation 
timescale r is increased, according to (9). (C,D) Bifurcation diagram showing supercritical 
{^-h) and subcritical (I^u) Hopf bifurcations that arise as the input is increased for (C) 
T = 10 and (D) t = 100. Firing rate function is sigmoidal (2). Other parameters are 0 = 1, 
a = 0.5, and 7 = 15 


Roots of (7), defining fixed points of (1) are plotted as a function of the input I 
in Fig. 1C,D. Utilizing the sigmoidal firing rate function / given by (2), we can 
show that there will be a single fixed point as long as <0 > a. In this case, we can 
compute 


dg(M) 

du 


-{(j) - a)f'{{a - (j))u + I) 


(<0 - a)e“T'd“-'^)®+^) 

.41_ i ____ ^ Q, 

_|_ e-7((a-0)ti+l))^ 


Since u is monotone increasing, then u — g{u) is monotone increasing. Further, 
noting lima-i-zboo [m — g{u)] = ± 00 , it is clear u — g{u) crosses zero once, so (7) has 
a single root when (p > a. Stability of this equilibrium is given by the eigenvalues 
of the associated Jacobian 


J{u, a) 


1 + «/'{(« - 4>)u + I) -/'((a - 4')u + I) 

(0/t “1/'^ 
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We note that the sigmoid (2) satisfies the Ricatti equation f' = 7/(1 — /), so we 
can use (7) to write 

T,- f—I + a^u(l — u) —ju{l — u)\ 

<A/r -1/r )■ 

Oscillations arise when stable spiral equilibria destabilize through a Hopf bifurca¬ 
tion. Hopf bifurcations will occur when complex eigenvalues associated with fixed 
points (m, a) cross from the left to the right half plane. We can require this with the 
pair of expression: tr(J) = 0 and tr(J)^ < 4det(J). Thus, a necessary condition 
for the Hopf bifurcation point is that the equilibrium value u satisfy 

a'yu{\ — u) = 1 -|- 1 /t. 


Solving this for u yields 

u±H = ^ [1 ± y/1 - 4x] , (8) 

Thus, Hopf bifurcations will only occur when the timescale of adaptation is suffi¬ 
ciently large t > [ay/d — 1]~^. Plugging the formula ( 8 ) back into the fixed point 
equation (7) and solving for the input I, we can parameterize Hopf bifurcation 
curves based upon the equation 


I±H = - In 

7 


U±H 

1 - U±H 


(a — (j))u±H, 


along with the additional condition tr(J)^ < 4det(J) which becomes 


( 9 ) 


4 4(/ A(j) 

ar^ ar 


( 10 ) 


which will always hold as long as </> > a. We partition the parameter space (r, I) 
using our formula for the Hopf curve (9) in Fig. IB. As demonstrated, there tend 
to be either two or zero Hopf bifurcation points for a given timescale t, and the 
coalescence of the two Hopf points is given by the point where t = [ay/d — 1]”^. 

In the limit of slow adaptation t ^ 1, we can separate the timescales of the 
activity u and adaptation a variables, finding u will equilibrate according to the 
equation 


u{t) = f{au{t) - a{t) + I), (11) 

and subsequently a will slowly evolve according to the equation 

d{t) = [4>u{t) — o] /r. ( 12 ) 

We always have an implicit formula for u{t) in terms of a{t), so the dynamics 
will tend to slowly evolve along the direction of the a variable. This demonstrates 
why periodic solutions to ( 1 ) are comprised of a slow rise and decay phase of 
a, punctuated by fast excursions in the activity variable u. In general, it is not 
straightforward to analytically treat the pair of equations ( 11 ) and ( 12 ), but we 
will show how computing solutions of the singular system becomes straightforward 
when we take the high gain limit j —)■ 00 . 
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Fig. 2 Analytical approximations to periodic solutions of (1) with a Heaviside firing rate func¬ 
tion (3). (A) Numerical simulation (solid lines) of the periodic solution is well approximated 
by the analytical approximation (dashed lines) given by (15) when I = 0.2 and r = 100. (B) 
The period of the oscillation T computed from numerical simulations (dots) is accurately ap¬ 
proximated by the analytical formula (solid lines) given by (14). Other parameters are a = 0.5 
and (p = 1. 


Having established the existence of periodic solutions to (1) in the case of 
sigmoid bring rates ( 2 ), we now explore the system in the high gain limit 7 —t oo 
whereby the bring rate function becomes a Heaviside (3). In this case, hxed points 
(u, a) satisfy the equations 


u = H{{a — 4>)u -h I) 


1 ■. u < I/{(j) — a), 
0 : u > I/{(/) — a), 


and a = cpu. Thus, assuming (j) > a, then u = 0 when / < 0 and u = 1 when I > 
{(j} — a). In both cases, the hxed points are linearly stable. When 0 < I < (0 — a), 
there are no fixed points and we expect to hnd oscillatory solutions. Assuming 
r ^ I, we can exploit a separation of timescales to identify the shape and period 
of these limit cycles. To begin, we note that on fast timescales 


u{t) = —u{t) + H{I + au{t) — oo), 

where ao is a quasi steady state. On slow timescales on the order of r, then u(t) 
quickly equilibrates and 


u{t) = H{I + au{t) — a{t)), (13a) 

Ta{t) = —a{t) + (13b) 

Periodic solutions to (1) must obey the condition (w(t),o(t)) = {u{t + nT),a(t + 
nT)) for t G [0,T] and n G Z, so we focus on the domain t G [0,T]. Examining 
(13), we can see oscillations in (1) involve switches between u(t) ss 1 and u{t) « 0. 
We translate time so that u{t) ss 1 on f G [0,Ti) and u{t) ss 0 on t G [Ti,T). 
Subsequently, this means for t G [0,Ti] the system (13) becomes u = 1 and 
rd = —a + 0 so a{t) = 0 — (0 — . We know a(0) = I because m(0“) = 0 

in (13), and the argument of H{x) must have crossed zero at f = 0. In a similar 
way, we hnd on t G [Ti,T) that u = 0 and a{t) = (/ + Using the 














Population models of up and down states 


9 


conditions a(Ti) = Ia and a(T) = I, we find that the rise time of the adaptation 
variable (or the duration of the up state) is 


Ti = T In 


i-/ 


< — a — I \ 

and the decay time (or the duration of the down state) is 

I + al 


T 2 = T In 


and the total period of the oscillation is 


T = T In 


(/ + a)(((>-/) 


(14) 


Thus, approximate periodic solutions to (1) in the case of a Heaviside firing rate 
(3) take the form 


r 1 e [0,Ti) 

\0 :te [Ti,T] 

r<^-(0 -7)e-‘/r :tG[0,Ti), 

1 (/ + a)e-(‘-^i)/" :te[Ti,T]. 


(15a) 

(15b) 


We demonstrate the accuracy of the approximation (15) in Fig. 2A. Furthermore, 
we show that relationship between the period T and model parameters is well 
captured by the formula (14). Notice there is a non-monotonic relationship between 
the period T and the input I. We can understand this further by noting that the 
rise time T\ of the adaptation variable a increases monotonically with input 


dTi 

dl 


ra 

{(j) - I){(j) -a- I) 


> 0 , 


when 0 < / < {(j) — a). Furthermore, the decay time T 2 of the adaptation variable 
a decreases monotonically with input 


dTa 

dl 


TO. 

I{I + a) 


< 0 , 


when 0 < I < (cj) — a). Thus, as / —> 0^, the slow oscillation’s period T is 
dominated by very long decay times T 2 ^ 1 and as /—)■((() — a) “, it is dominated 
by very long rise times Ti ^ 1. We can identify the minimal period as a function 
of the input I by finding the critical point of T{I). To do so, we differentiate and 
simplify 


dT Ta4>{21 — {(f> — a)) 

Tf ^ ~ I{I + - a - I)’ 

so the critical point of T(7) is Icrit = (<('~tt)/2, which corresponds to the minimal 
value of the period Tmin{I) = 2Tln [{(j) -)- a)/{(j) — a)] as pictured in Fig. 2B. 
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Fig. 3 (A, B, C) Periodic solution {u,a) and (D, E, F) phase sensitivity function 
of (1) plotted as a function of phase 0 = t/T for a sigmoidal firing rate function (2). (A,D) 
For shorter adaptation timescale r = 10 and input I = 0.2, the activity variable u has a more 
rounded trajectory, so perturbations to activity influence the oscillation phase more (note size 
of lobes in on Zu in (D). (B,E) As the adaptation timescale is increased to r = 100, with 
I = 0.2, the influence of perturbations to the activity variable decrease (compare lobes of 
to those in (D). Perturbations of the adaptation variable influence the phase more strongly 
as shown by the change in the relative amplitude of Za- (C,F) Increasing the input I = 0.4, 
with T = 10, increases the relative duration of the rise time of a. As a result, there is a wider 
region where perturbations to a advance the phase. Other parameters are a = 0.5, 0 = 1, and 
7 = 15. 


4 Phase response curves 

We can further understand the dynamics of the slow oscillations in (1) by compnt- 
ing phase response curves for both the case of a sigmoidal firing rate ( 2 ) and the 
Heaviside firing rate (3). As we will show, perturbations of the activity variable u 
have decreasing impact as the timescale of adaptation r and the gain 7 of the firing 
rate are increased. Pertnrbations of the adaptation variable o tend to dominate 
the resulting dynamics, as it is the evolution of this slow variable that primarily 
determines the phase of the oscillation. 

To begin, we derive a general formula that linearly approximates the influence 
of small perturbations on limit cycle solutions {uo{t),ao{t)) to (1). Essentially, we 
utilize the fact that solutions Z(t) to the adjoint equation associated with lineariza¬ 
tion about the limit cycle solution (uo(t), ao{t)) provide a complete description of 
how infinitesimal perturbations of the limit cycle impact its phase (Ermentrout, 
1996; Brown et ah, 2004). To start, we note that 

„ / ui \ _ / ui + ui - af'{auQ - uq + I)ui -b f'{auQ - uq + /)oi,\ _ /0\ 

V^aiy l^a'i - -b ai/r J V^/ 

is the linearization of (1) about the limit cycle (uo{t),ao{t)). Defining the inner 
product on T-periodic functions in as {F{t),G{t)) = F{t) ■ G{t)dt, we can 
find the adjoint operator £.* by noting it satisfies {F,CG) = (jC*F,G) for all 
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integrable vector functions F, G. We can then compute 

_ / -i) + v - af'(auo - ao + I)v - (^hjr 
\b J \-b+f'{auo-ao +I)v + b/T 

It can be shown that the nnll space of C* describes the response of the phase of 
the limit cycle {uo{t),ao{t)) to infinitesimal perturbations (Brown et ah, 2004). 
Note that if (uo(t), ao(t)) is a stable limit cycle then the nullspace of £ is spanned 
by scalar multiples of (uQ(f), aQ(f)). Furthermore, appropriate normalization re¬ 
quires that Z(t) ■ (uo(t), ao(t)) = 1 along with C*Z = 0 (Ermentrout, 1996). To 
numerically compute Z(t) = {Zu{t), Za{t)), we thus integrate the system 

Zu{t) = Zu{t) - af'{auo{t) - ao{t) + I)Zu(t) - (pZaifij/r, (17a) 

Za{t) = f'{auo{t) - ao(t) + I)Zu{t) - Za{t)fT, (17b) 



backward in time, taking the long time limit to find {Zu{t), Za{t)) on t G [0,T], 
and normalizing {{Zu{t), Za{t)), {uo(t), aQ{t))) = 1 by rescaling appropriately. We 
demonstrate this result in Fig. 3, showing the relationship between the shape and 
relative amplitude of the phase sensitivity functions {Zu, Za) and the parameters. 
Notably, perturbing the activity variable u become less influential as the timescale 
of adaptation t is increased {Zu)- Furthermore, there is a sharper transition be¬ 
tween phase advance and phase delay region of the adaptation phase response {Za) 
for larger timescales r. 

In addition to a general formula for the phase sensitivity functions {Zu{t), Za{t)), 
we can derive an amplitude dependent formula for the response of limit cycle so¬ 
lutions {uo{t),ao{t)) of (1) with a Heaviside firing rate (3), assuming t ^ 1. In 
this case, we utilize the formula for the period (14) and limit cycle (15), derived 
using a separation of timescales assumption. Then, we can compute the change 
to the variables (u, a) as a result of a perturbation (6u,ba), which we denote 

{uo{t), ao{t)) {uo{t),ao{t)). We are primarily interested in how the relative 

time in the limit cycle is altered by a perturbation 6u - how much closer or further 
the limit cycle is to the end of the period T after being perturbed. We can readily 
determine this by hrst inverting the formula we have for (Mo(t), oo(t)), given by 
(15), to see how this value determines the time to along the limit cycle 


to (mo, ao) 


Tin [{(/)-I)/{c/> - ao)] : MO = 1, 

Tin [{(j) - I){I + a)/ao/{(j) - a - /)] : mq = 0. 


(18) 


Using this formula, we can now map the value (mo,5o) to an associated updated 
relative time to along the oscillation. 

Here, we decompose the impact of perturbations to the u and a variables. We 
begin by studying the impact of perturbations 5u to the activity variable u. We 
can directly compute 


uo{t) = H{I + a [uo{t) + 5u] - ao{t)). 

Thus, the singular system (13) will be unaffected by such perturbations if sgn(/ + 
a [m -h 5u] —a) = sgn(/ -|- au — a) . This is related to the flatness of the susceptibility 
function Zu over much of the time domain in Fig. 3D,E,F. However, if sgn(/ + 
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Fig. 4 Phase response curves of the fast-slow timescale separated system r 1. (A,B) 
Amplitude Su- and 5a-dependent phase response curves Gu(^, ^u) and Ga(^, ^u) characterizing 
phase advances/delays resulting from perturbation of neural activity u and adaptation a. We 
compare analytical formulae (solid lines) to numerically computed PRCs (dashed lines). (C) 
Phase response curve associated with perturbations of the adaptation variable a in the small 
amplitude 0 < |5a| "C 1 limit. We compare the large amplitude formula (solid line) determined 
by (24) to the linear approximation (dotted line) given by (25) to numerical computations 
(dashed line). 


a[u + Su] — a) 7 ^ sgn{7 + aM —a), then uo{t) = l — uo{t), as detailed in the following 
piecewise smooth map: 


uo{t) 

uo{t) 

uo(t) 

uo{t) 


0 !-)■ uo(i) = 1 
0 1 -^ uo{t) = 0 

1 !->■ UQ(t) = 0 

1 !->• uo(t) = 1 


: Su > -{I - ao{t))/a > 0, 

: Su < -{I - ao{t))/a > 0 , 

: -Su < -{I + a - ao{t))/a < 0 , 
: -Su > -{I + a - ao{t))/a < 0 , 


where {uo{t),ao{t)) are defined by (15). The formnla (18) can then be ntilized to 
compnte the updated relative time to := lo(wo, ho); finding 

{ Tin [{(j) - /)(/ + a)/ao/{4i - a - I)] : Su > {I + a - ao)/a >0, mq = 1 
T + Tin [((^ - /)/(<^ - ao)] : -Su > (ao - 1)1 a >0, uo = 0 

to(uo,ao) : otherwise, 

(19) 

where ao = 4> — {<P — if uq = 1 and ao = (/ + if uo = 0. We 

can refer to the function tojT, where to is defined by (19), as the phase transition 
curve for u perturbations. Thus, the function Gu{d,Su) = (to ~ to)/T will be 
the phase response curve, where 0 = tojT, and phase advances occur for positive 
values and phase delays occur for negative values. We plot the function Gu(0,Su) 
in Fig. 4A for different values of Su, demonstrating the nontrivial dependence 
on the perturbation amplitude is not simply a rescaling but an expansion of the 
non-zero phase shift region. Due to the singular nature of the fast-slow limit cycle 
(15), the size of the phase perturbation has a piecewise constant dependence on the 
amplitude of the u perturbation. Note, this formulation allows us to quantify phase 
shifts that would not be captured by a perturbative theory for phase sensitivity 
functions, as computed for the general system in (17). 

For perturbations Sa of the adaptation variable a, there is a more graded de¬ 
pendence of the phase advance/delay amplitude on the perturbation amplitude Sa- 
We expect this, as it was a property we observed in Za as we varied parameters in 
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Fig. 3. We can partition the limit cycle {uo{t), ao{t)) into fonr different regions: two 
advance/delay regions of exponential satnration and two early threshold crossings. 
First, note if Uo{t) = 1 and ao{t) + 5a < / + a, then 

Uo(t) = 1, aQ{t) = 4) - {(j) - + 5a, (20) 

so to = Ti-tyj with tu, = Tin [{4> - no - 5a)/{(l) - I - «)], bnt if ao(t) + 5a > I + a, 

then 

uo(t) = 0, ao{t) = (j) - {(j) - + 5a- (21) 

Determining the relative time of the perturbed variables (uo(t), oo(t)) in (20) is 
straightforward using the mapping (18). However, to determine the relative time 
described by (21), we compute the time, after the perturbation, until do{t) = I+a, 
which will be tiu = TIn [(no + 5a)/{I + a)], so to = Ti — tw Second, note if uo{t) = 
0 and ao(t) + 5a > I, then 

uo(t) = 0, ao(t) = (/ + + 5a, (22) 


so to = T — tw with to] = Tin [(no + 5q)/I], but if oo(t) + 5a < I, so that it is 
necessary that 5a < 0 , then 

uo(t) = 1 , do(t) = (I + + Sa- (23) 

In the case of (23), we note that to, = T\n[{(j)— oq — 5a)/{4>— I)], so to = T — tw 
Combining our results, we hnd we can map the relative time to the perturbed 
relative time as 


' Ti - T In [(0 - ao - 5a)/(()) - / - a)] : 5a < I + a - ao, uq = 1, 

Fi - Tin [(ao + 5a)/(/+ a)] : 5a > / + a - oo, uo = 1, 

to = < (24) 

T - Tin [(no + 5a)//] : 5a > / - no, uo = 0, 

_ T - T In [(0 - ao - 5a)/ (</-/)] : 5a < I - ao, uo = 0, 

where ao = (/—(</ — I)e~*°^'^ if uo = 1 and ao = (/ + if uo = 0 . 

Again, we have a phase transition curve given by the function io/T and phase 
response curve given by Ga{0,5a) = {to — to)/T, where 6 = to/T. As opposed to 
the case of u perturbations, the phase perturbation here depends smoothly on the 
amplitude of the a perturbation 5a- 

Furthermore, we can obtain a perturbative description of the phase response 
curve for the singular system (13) in two ways: (a) Taylor expand the amplitude- 
dependent phase response curve expressions dehned by (19) and (24) and truncate 
to linear order or (b) solving the adjoint equation (17) in the case of a Heaviside 
hring rate (3) and long adaptation timescale t ^ 1- We begin with the first 
derivation, which simply requires differentiating (19) to demonstrate that the in¬ 
finitesimal phase response curve (iPRC) associated with perturbations of the u 
variable is zero almost everywhere. However, differentiating (24) reveals that the 
iPRC associated with perturbations of the adaptation variable a is given by the 
piecewise smooth function 


tjr 

T{4>-I)s'' 

+t-T+/r 

T{I + a) 


■-1 e (o,Ti), 

e (Ti,T). 


( 25 ) 
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Furthermore, note we could derive the same result by solving the adjoint equations 
(17) in the case of Heaviside firing rate (3), so that 

Zu{t) = -Zu{t) + a5{auo(t) - ao{t) + I)Zu{t) + (j)Za{t)/T, (26a) 

Za{t) = -6{auo{t) - ao{t) + I)Zu{t) + Za{t)fT. (26b) 

Note, we have reversed time t i— —t, so we can simply solve the system forward. 
Furthermore, we can use the identity 

Utilizing the separation of timescales, t ^ 1, we find that almost everywhere 
(except where t = 0,Ti,T), we have that (26) becomes the system 

Zu{t) = -Zu{t), TZa{t) = Za{t). (28) 

As before Zu{t) will be zero almost everywhere, whereas Za(t) = , where 

A{t) is a piecewise constant function taking two different values on t G (0, Ti) and 
t G (Ti,T), determined by considering the 6 distribution terms. This indicates 
how one would derive the formula (25) using the adjoint equations (26). 

Note, in previous work (Jayasuriya and Kilpatrick, 2012), we explored the en¬ 
trainment of slowly adapting populations to external forcing, comprised of smooth 
and non-smooth inputs to the system (1). In the next section, we explore the im¬ 
pact of external noise forcing on the slow oscillations of (1), subsequently demon¬ 
strating that noise can be utilized to entrain the up and down states of two distinct 
networks. 


5 Impact of noise on the timing of up/down states 

We now study the effects of noise on the duration of up and down states of the 
single population model (1). Switches between high and low firing rate states 
can occur at irregular intervals (Sanchez-Vives and McCormick, 2000), suggesting 
internal or external sources of noise determine state changes. This section focuses 
on how noise can reshape the mean duration of up and down residence times. Due 
to our findings in the previous sections, we focus on noise applied to the adaptation 
variable in this section. As we have shown, very weak perturbations to the neural 
activity variable have a negligible effect on the phase of oscillations. Analytic 
calculations are presented for the piecewise smooth system with Heaviside firing 
rate (3), as accurate approximations of the mean up and down state durations can 
be computed. 

Our approach is to derive expressions for the mean first passage times of both 
the up and down state (Ti and T 2 ) of the stochastic population model (4). Focusing 
on adaptation noise allows us to utilize the separation of fast-slow timescales, and 
recast the pair of equations as a stochastic-hybrid system 

u(t) = H({au{t) + I — a{t)), 
da{t) = [-a{t) -h dt/r -|- dfa(f). 
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Fig. 5 Noise alters the duration of up and down states. (A) Numerical simulation of the 
stochastically driven population model (4) demonstrates up and down state durations (e.g., 
T\ and T 2 ) are variable when driven by adaption noise with = cr^f, Ua = 0.01. Switches 
are determined by the threshold crossings of the adaptation variable a{t) = I and a{t) = 7 + 0 :. 
(B) Up/down states become more variable when the noise amplitude Ga = 0.02. (C) Mean 
durations of the up and down state, (Ti) and {T 2 ), decrease as a function noise amplitude 
cTfl. (D) Impact of noise Ga on the balance of up to down state durations T 1 IT 2 as input I 
is varied. Firing rate is given by the Heaviside function (3). Other parameters are o = 0.5, 
(p = 1, and T = 50. 


where is white noise with mean {^a) = 0 and variance To begin, 

assume the system has just switched to the up state, so the initial conditions are 
u{0) = 1 and a(0) = /. Determining the amount of time until a switch to the down 
state requires we calculate the time Ti until the threshold crossing a(Ti) = I a 
where a{t) is determined by the stochastic differential equation (SDE) 

do(t) = [-a{t) + (j)] dt/r + d^a, 

which is the well-known threshold crossing problem for an Ornstein-Uhlenbeck 
process Gardiner (2004). The mean Ti of the passage time distribution is thus 
given by defining the potential V{a) = ^ ^ and computing the integral 

1 /*/+a rx 

Ti = 4 

JI J — oo 

Next, note that the duration of the down state T 2 will be the amount of time until 
the threshold crossing o(T 2 ) = I given u(0) = 0 and a(0) = I + a, where a(t) 
obeys the SDE 


do(t) = [-a(t)] dt/r -|- d§a(t)- 
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2 

Again, defining the potential (a) = , we can compute 

J — I —a J— oo 

We compare the theory we have developed utilizing passage time problems to 
residence times computed numerically in Fig. 5C. Notice that increasing the noise 
amplitude tends to shorten both up and down state durations on average, due to 
early threshold crossings of the variable a{t). 

Furthermore, we can examine how noise reshapes the relative balance of up 
versus down state durations. Specifically, we will explore how the relative fraction 
of time the up state persists T'i/(T’i + T 2 ) changes with noise intensity aa and 
input I. First, notice that, in the absence of noise the ratio 

Ti ^ ln[(0-7)/((A-n-/)] 

T 1 +T 2 ln[(/ + a)(<^-7)//(<^-a-/)]' ^ ’ 

The up and down state have equal duration when Ti/(Ti + T 2 ) = 1/2, or when 
the input I = {<p — a)/2, as shown in Fig. 5D. Interestingly, this is the precise 
input value at which the period obtains a minimum, as we demonstrated in section 
3. Along with our plot of (29) in the noise-free case (era = 0), we also study the 
impact of noise on this measure of up-down state balance. Noise leads to up and 
down state durations becoming more similar, so the ratio (29) of the means Ti and 
T 2 flattens as a function of the input 7. This is due to the fact that long durations, 
wherein the variable a{t) occupies the tail of exponentially saturating functions 
Aq + , are shortened by early threshold crossings due to the external noise 

forcing. 


6 Synchronizing two uncoupled populations 

Now we demonstrate that common noise can synchronize the up and down states 
of two distinct and uncoupled populations. We begin with the case of identical 
noise and then, in section 7, relax these assumptions to show that some level of 
coherence is still possible when each population has an intrinsic and independent 
source of noise. This is motivated by the observation that the neural Langevin 
equation derived in the large system-size limit of a neural master equation tends 
to possess intrinsic noise in each population, in addition to an extrinsic common 
noise term (Bresslofl and Lai, 2011). As we will show, intrinsic noise tends to 
disrupt the phase synchronization due to extrinsic noise. 

To begin, we recast the stochastic system (5), describing a pair of adapting 
noise-driven neural populations, as a pair of phase equations: 

d0i(t) = wdt-h Z(6i(t)) • d^(t), (30a) 

d6>2(t) = a;dt-hZ(6»2(t)) •d^(t), (30b) 

where 61 and 02 are the phase of the first and second neural populations. As 
we demonstrate in Fig. 6A, this introduction of common noise tends to drive 
the oscillation phases 0i{t) and 02 {t) toward one another. Note that since the 
governing equations of both populations are the same, then the phase sensitivity 
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Fig. 6 Synchronizing slow oscillations in two uncoupled populations described by (5) with 
sigmoidal firing rate (2). (A) Single realization of the system (5) driven by common noise 
to the adaptation variable ((^q) = e'^t, e = 0.01) with input I = 0.2 and adaptation timescale 
T = 50. Notice that the phase difference ip{t) = Ai{t) — A 2 {t) roughly decreases over time. 
(B) Plot of the log of the phase difference y(t) = for several realizations (thin lines) 

compared with the theory (thick line) of the mean y(0) + \t computed using the Lyapunov 
exponent (35). (C) Lyapunov exponent A decreases as a function of the adaptation timescale 
r, for I = 0.2. We compare numerical simulations (dots) to theory (solid). (D) Lyapunov 
exponent A varies non-monotonically with the strength of the input I. Other parameters are 
a = 0.5, 7 = 15, and (f) = 1. 


function Z(^) will be the same for both. Furthermore, the synchronized solution 
Oi{t) = 02 {t) is absorbing - once the phases synchronize, they remain so. We 
can analytically calculate the Lyapunov exponent A of the synchronized state to 
determine its stability. In particular, we will be interested in how this stability 
depends on the parameters that shape the dynamics of adaptation. 

We next convert the pair of Stratonovich differential equations into a equivalent 
pair of Ito differential equations: 


dOiit) = [w + Z'(6'i(t))'^DZ(6»i(i))] dt + ■ d^{t), (31a) 

d02{t) = [w + Z'(02{t))'^DZ(6»2(t))] dt + Z{02{t)) ■ d${t), (31b) 

introducing a drift term due to our change in definition of the noise term. Now, 
to determine stability of the synchronized state Oi{t) = 02 (t), we assume we are 
inhnitesimally close to it. We dehne tpit) = ^i(f) ~ ^ 2 ( 1 ) and require \ip(t)\ ^ 1- 
Linearizing the system of Ito differential equations (31) with respect to ipit) then 
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yields 

d^{t) = ^{t) (z'{9{t)f-DZ{e{t))ydt + z'{e(t))-d^ , (32) 

where 0{t) obeys either one of the equations in (31). Applying the change of 
variables y{t) = In-0(1), we can rewrite (32) as 

dy{t) = (Z'DZ)' dl - (^Z'^DZ'j dl + Z' • d^(l). (33) 

Notice, on average, the log of the phase difference y{t) tends to decrease over time 
(Fig. 6B), indicating the phases 6 i and 02 move toward one another. Subsequently, 
we can integrate equation (33) to determine the mean drift of y{t) 

r t 

A:=^lhn j [(Z'(6»(s))DZ(6>(s)))'- (^Z''^(6»(s))DZ'(6>(s))^] ds, 

so the phase difference 0(1) will tend to decay grow if the Lyapunov exponent 
A < 0 (A > 0), and the synchronous state will be stable (unstable). Now, utilizing 
ergodicity of (33), we can compute A using the ensemble average across realizations 
of Z'(0{t)) ■ ^(1), so 

\ = I Ps{0) [(Z'(0)DZ(d))' - (z'^(e)DZ'(e))] d9, (34) 

where Ps{0) is the steady state distribution of 0. Since noise is weak (Djfc ^ 1, 
j,k = 1,2), we can approximate the distribution as constant Ps{0) = 1. Upon 
applying this to the integrand of (34) and noting the periodicity of Z{0), we hnd 
we can approximate the Lyapunov exponent 

A = - / Z'^(0)DZ'(0)d0. (35) 

./o 

Assuming noise to the activity variable u and adaptation variable a is not cor¬ 
related, D will be diagonal. In this case, we can further decompose the phase 
sensitivity function into its Fourier expansion 

OO 

m = E Rfc sin(27rfc0) -|- cos(27rfc0), 
k=0 

where = (afci,afc 2 )^ and = (bfci,bfc 2 )^ are vectors in so that 

OO 

z'(0) = E 27rfe [afc cos(27rfc0) — b^ sin(27rfc0)] , 
fc =0 

and we can expand the terms in (35) to yield 

OO 

A = — + bfcij Dn -|- {^^2 + ^k2'j D22^ ■ 

k=0 

Thus, as long as Z{0) is continuous and non-constant, the Lyapunov exponent A 
will be negative, so the synchronous state 0i = 02 will be stable. Note, continuity is 
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Fig. 7 Stationary density Mo{ip) of the phase difference ■i/i = 6 i — 62 for two slowly oscillating 
neural population driven by both common and independent noise ( 6 ). As the degree of noise 
correlation is decreased from (A) Xa = 0.95 to (B) Xa = 0.90, the density spreads, but there 
is still a peak at i/j = 0, the phase-locked state. We focus on noise in the adaptation variable, 
so (Tu = 0 and era = 0.01. Other parameters are a = 0.5, 7 = 15, 0 = 1, and r = 20. 


not satisfied in the case of our singular approximation to Z(0). We demonstrate the 
accuracy of our theory (35) in Fig. 6C,D, showing that A decreases as a function 
of T and is non-monotonic in I. Thus, slow oscillations with longer periods are 
synchronized more quickly, relative to the number of oscillation cycles. Since the 
Lyapunov exponent has highest amplitude |A| for both low and high values of the 
tonic input I, we also suspect this is related to the period of the oscillation T. 


7 Impact of intrinsic noise on stochastic synchronization 

We now extend our results from the previous section by studying the impact 
of independent noise in each population. Independent noise is incorporated into 
the modified model ( 6 ). Noting, again there is a periodic solution to the noise-free 
version of this system, phase-reduction methods can be used to obtain approximate 
Langevin equations for the phase variables (Nakao et ah, 2007) 

d6i = uidt + Z(0i(t)) • [d$c(^) + d$i{f)] , (36a) 

d02 = uidt + Z(02(t)) • [d$c(^) + d$ 2 (i)] I (36b) 

where the noise vectors = {Xu^uc, Xa^ac)’^ and = {^/l- xl^uj, \/l- xHaj)'^ 
{j = 1,2). We can reformulate the pair of Stratonovich differential equations as 
Ito stochastic differential equations given by the system 

dOi = Ai{e)dt + dCi{e,t), (37a) 

d02 = A2{e)dt + dC2{0,t), (37b) 

where the statistics of the noise terms dCj{6,t) = Z( 0 j(t)) • [d^c(^) + U — 

1 , 2 ) are specified by {d(^j{6,t)} = 0 and {dCj{d,t)dCk{S,t)} = Cjk{0)dt where 

Cjk{0) = (xuTtii^Zuiej) + Xa-Dai^ZaiOj)) (^Xu-Dli^ZuiOk) + Xa-Dli^ZaiOk)) 


(38) 
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separating the impact of correlated and local sources of noise. The drift terms 
can thus be calculated Aj(0) = w + \ The Fokker-Planck equation 

describing the evolution of the probability density function P{0,t) of the phases 
is given 

f = - E + l±±^ ■ (39) 

j = l j = lk = l 

Now, we apply a change of variables to the Fokker-Planck equation (39) dehned 
6 j = Assuming noise is weak, the function <5(1?, t) varies slowly compared 

with the period of the phase oscillators 6j. Thus, we average the drifts Aj(0) 
and correlation function Cjk{0) over a single period. The resulting Fokker-Planck 
equation is then 


dQ{'0,t) 

dt 


1 ^ 2 

iEE 


a" 

ddjd'&k 


\Bjk{ 0 )Q] 


where the averaged correlation function is given by the formula 


Bjk{'0)=g{0i-e2) + h{Q)6j,k, 


where the correlation functions are dehned 

9{0) = £ [xl'DucZu{0')Zu{9' + 9) + xl'DacZa{9')Za{9' + 0 )] d9' 

and 

h{9) = £ [(1 - xl)'DuiZu{9')Zu{9' + 0) + (1 - xl)B>aiZa{9')Za{9' + 9)] d9'. 

We study the relationship between the phases i9i and 02 by substituting the for¬ 
mula for the averaged correlation matrix 


dQ{'d,t) 

Wt 


^ [ 5 ( 0 ) -h h( 0 )] 


'd'^Q d^Q' 
d 0 l ^ 30^ 


d 0 id 02 


[c/(r?i - 02 )Q] ■ 


(40) 


We can write (40) as a separable equation by employing a change of variables that 
tracks the average p = (iJi + 02)/2 and phase difference tp = 0i —02 oi the original 
position variables, so 


dU{p,t) 

dt 

dt 


[ 5 ( 0 ) - gW + h( 0 )] M(^, t). 


(41a) 

(41b) 


Thus, we can solve for the stationary solution of the system (41a) by serving Ut = 
Mt = 0 and requiring periodic boundary conditions. We hnd that the stationary 
distribution of the position average is Uo{p) = 1. In addition, we can integrate the 
stationary equation for M(%l),t) to hnd 

_mo_ 

(^1 [(2 - xS)fl«(0) - xlguii’)] + CTa [(2 - Xa)fla(O) - Xa5a(^)] ’ 


Mo(^) 


(42) 
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where mo = 1/ Jq [fl(0) — g(x) + h( 0 )] ^ da; is a normalization factor and we have 
simplified the expression using Gui = Du 2 = = cr^ and Dai = Da 2 = Da/ = 

<Ja and defined 

r z,{e)z,{e + ^)de. 

Jo 

When noise to each layer is independent {xu, Xa —>■ 0), then Moitp) = 1 is constant 
in space. When noise is totally correlated (xu^Xa —> 1), then Mq(% 1)) = 5((j)). 
The stationary distribution Moixp) will broaden as the correlations Xu and Xa 
are decreased from unity, with a peak remaining at 0 = 0. We demonstrate the 
accuracy of the formula (42) for the stationary density of the phase difference 
in Fig. 7, showing that it widens as the level noise correlation is decreased. Again, 
we focus on the impact of adaptation noise. Thus, even when independent noise 
is introduced, there is some semblance of synchronization in the phases of two 
noise-driven neural populations (6). 


8 Discussion 

We have studied the impact of deterministic and stochastic perturbations to a 
neural population model of slow oscillations. The model was comprised of a sin¬ 
gle recurrently coupled excitatory population with negative feedback from a slow 
adaptive current (Laing and Chow, 2002; Jayasuriya and Kilpatrick, 2012). By 
examining the phase sensitivity function {Zu, Za), we found that perturbations of 
the adaptation variable lead to much larger changes in oscillation phase than per¬ 
turbations of neural activity. Furthermore, this effect becomes more pronounced as 
the timescale t of adaptation is increased. Introducing noise in the model decreases 
the oscillation period and helps to balance the mean duration of the oscillation’s 
up and down states. When two uncoupled populations receive common noise, their 
oscillation phases 6 \ and O 2 eventually become synchronized, which can be shown 
by deriving a formula for the Lyapunov exponent of the absorbing state = O 2 
(Teramae and Tanaka, 2004). When independent noise is introduced to each popu¬ 
lation, in addition to common noise, the long-term state of the system is described 
by a probability density for tp = 9i — 62 , which peaks at tp = 0 . 

Our study was motivated by the observation that recurrent cortical networks 
can spontaneously generate stochastic oscillations between up and down states. 
Guided by previous work in spiking models (Compte et ah, 2003), we explored a 
rate model of a recurrent excitatory network with slow spike frequency adaptation. 
One of the open questions about up and down state transitions concerns the degree 
to which they are generated by noise or by more deterministic mechanisms, such 
as slow currents or short term plasticity (Cossart et ah, 2003). Here, we have 
provided some characteristic features that emerge as the level of noise responsible 
for transitions is increased. Similar questions have been probed in the context 
of models of perceptual rivalry (Moreno-Bote et ah, 2007). In addition, we have 
provided a plausible mechanism whereby the onset of up and down states could 
be synchronized in distinct networks (Volgushev et ah, 2006). 

There are several potential extensions to this work. For instance, we could ex¬ 
amine the impact of long-range connections between networks to see how these 
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interact with common and independent noise to shape the phase coherence of 
oscillations. Similar studies have been performed in spiking models by Ly and Er- 
mentrout (2009). Interestingly, shared noise can actually stabilize the anti-phase 
locked state in this case, even though it is unstable in the absence of noise. Fur¬ 
thermore, it is known that coupling spanning long distances can be subject to 
axonal delays. In spite of this, networks of distantly coupled clusters of cells can 
still sustain zero-lag synchronized states (Vicente et ah, 2008). Thus, we could also 
explore the impact of delayed coupling, determining how features of phase sensi¬ 
tivity function interact with delay to promote in-phase or anti-phase synchronized 
states. 
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